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Abstract 

The Dyson-Schwinger equations arising from minimizing the vacuum energy density in the Hamil- 
tonian approach to Yang-Mills theory in Coulomb gauge are solved numerically. A new solution 
is presented which gives rise to a strictly linearly rising static quark potential and whose existence 
was previously observed in the infrared analysis of the Dyson-Schwinger equations. For the new 
solution we also present the static quark potential and calculate the running coupling constant 
from the ghost-gluon vertex. 
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I. INTRODUCTION 



Recently, there has been renewed interest in the Hamiltonian approach to Yang-Mills theory 
in Coulomb gauge [l|. In this gauge the Yang-Mills Schrodinger equation was solved varia- 
tionally in the continuum using Gaussian types of wave functionals 0, 0, 0] . Minimization 
of the energy density results in a set of Dyson-Schwinger equations. These equations have 
been solved analytically in the infrared using power law ansatze for the corresponding Green 
functions {5J. Thereby, two infrared solutions of the Dyson-Schwinger equations have been 
found which differ in the infrared exponents. Only one of these two solutions was previously 
found in the numerical solution of the Dyson-Schwinger equations in the whole momentum 
regime [3|. We present here the other one which has the attractive property that it gives 
rise to a strictly linear rising static quark potential. 

The organization of the paper is as follows: In the next section we briefly summarize the 
basic ingredients of the Hamiltonian approach to Yang-Mills theory in Coulomb gauge and 
present the resulting Dyson-Schwinger equations. The numerical solutions of these equations 
are presented in section 3. We also calculate from these solutions the static quark potential 
and the running coupling constant. Finally, our conclusions are given in section 4. 



II. THE DYSON-SCHWINGER EQUATIONS IN THE HAMILTONIAN AP- 
PROACH TO YANG-MILLS THEORY IN COULOMB GAUGE 



In Coulomb gauge, where d ■ A = 0, the Yang-Mills Hamiltonian reads pj 

H=\j J-'TUTl + ~ j 5 2 + £ J J-'p{-bd)-\-d 2 ){-bd)- l Jp . (1) 

Here, n"(x) = 5/i5A^(x) denotes the momentum operator, representing the color electric 
field, B is the color magnetic field and p a (x) = —A° b {x)T\ b {x) is the non-Abelian color charge 
of the gauge field. Furthermore, J(A) = Det(-Didi) is the Faddeev-Popov determinant with 
jjab _ fiabg. _j_ j^ab k e j n g the covariant derivative in the adjoint representation of the gauge 
group (A1 b = f acb Al with the structure constants f abc ). The first term in the Hamiltonian is 
the Laplacian in the curved space of gauge orbits, the second term represents the potential 
and the last term arises from the longitudinal momentum part of the kinetic energy after 
resolving Gauss' law. We use the following ansatz for the vacuum wave functional [31] 



nA] =N' r (-\jA«A) 



(2) 



where oj(\x — x'\) is a variational kernel determined by minimizing the energy 

(V\H\V) = J VAJ[A}^*[A]H^[A] . (3) 
This yields a set of coupled Dyson-Schwinger equations for the ghost propagator 

G ab (x,x') = (^Kx^K-Dd)- 1 ^'^)^) = -{-d^dix.x') , (4) 
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the gluon propagator 



D$M = (9\AS(x)A>&)\*) = l -8 a %{x)u-\x,x') , t^x) = % - ^ 
and the Coulomb form factor / defined by 

(^K-Ddy'i-d^i-Dd)- 1 ^} = {^(-Ddy^i-d^f^K-Dd)- 1 ^} 

A further quantity arising in the evaluation of (^\H\^) is the scalar curvature \i 



X(x,x r ) 



2{N? 



'tij(x)Xij {%■, X ) , Xij i x t x ) 



c 



6 2 In J [A] 



5A?(x)8A b Ax') 



(5) 



(6) 



(7) 



which denotes the ghost loop contribution to the gluon energy. This quantity measures the 
"curvature" of the space of transversal gauge orbits and vanishes in QED. 
Following ref. [jjj], we will solve here the Dyson-Schwinger equations to one-loop order. To 
this order it is sufficient to calculate the Coulomb form factor / perturbatively to ensure 
the proper ultraviolet asymptotics of the corresponding loop integrals (see ref. [3| for more 
details). In the infrared, this results in f(k — > 0) = const. After regularization and renor- 
malization using subtracted equations with some arbitrary renormalization scale fi, the re- 
maining Dyson-Schwinger equations for the ghost and gluon propagators and the expression 
for the curvature in momentum space read [3j 



d-\k) = d~ l {n) - AI d (k,/i) 



(8) 



cu 2 (k) 



X 



\k) + k 2 + 2 x (k) (AI« (*, fj) - A/« (0, n) + c) + A/^) (k, /x) + c , (9) 



(10) 



with c = cu 2 (/i) — X 2 (/ i ) ~ V? ~ 2x(/i)(— Alw(0) + c). Here <i _1 ( / u), u)(fi), x(/f)> an d c are 
renormalization constants. Furthermore, the AI X refer to the following differences of loop 
integrals 



AI d {k,/i) 
Al x (k,/i) 



N 

I d (k) - I d (n) , I d (k) = 

i x (k) - im , i x (k) = ^ 

4 n \k)-lM(ji),n = l,2 



d 3 q 
(2^)3 
d 3 q 
(2^ 



-"(k-q) 2 
l-(k-q) 2 



d(|k-q|) 
(k - q) 2 u(q) 
d(\k-q\)d(q) 
(k - q)2 



(11) 
(12) 



Nc 
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d 3 q 
(2tt)s 



l + (k-q) 



d(\k - q|) 2 /(|k - q|) Kg) - X (q)} n - M*) - X (*)] n 



(k-q) 



u;(g) 



(13) 



For the solution of the coupled Dyson-Schwinger equations (jSJ) and (Q, we implement the 
horizon condition 0] 

tT 1 ^ = 0) = . (14) 
Assuming power law behaviour for the Green functions in the infrared 



w(Jfe) 



, d(k) 



B 



(15) 
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the infrared analysis of the ghost Dyson-Schwinger equation yields the sum rule 

2p = a + l, (16) 

due to the non-renormalization of the ghost-gluon vertex, a feature of both Landau 0] 
and Coulomb gauge jl, H|. The horizon condition (JTJJ implies (3 > 0. For (3 > 1/2 we 
obviously obtain an infrared divergent gluon energy which is a manifestation of confinement. 
Furthermore, /3 > 1/2 also implies an infrared divergent curvature x(^)- Indeed, an infrared 
analysis jl] of Eq. (1T2"1) shows that for k — > 0, x(A;) oc k 213 " 1 which by the sum rule (lT6"j) is 
the infrared behaviour of u(k). This is consistent with the gap equation ([9]). For an infrared 
divergent x(k) the gap equation Q reduces in the infrared to 

lim(u(k)- X (k)) = c, (17) 

implying that u(k) and x{k) have the same infrared (divergent) behaviour. Eq. (iT7j) along 
with the ghost Dyson-Schwinger equation (jSJ) in the infrared limit can be solved analytically 
for the infrared exponents. It was found earlier that, using angular approximation, one 
finds only one solution f3 = 1 0] . However, performing an infrared analysis without angular 
approximation, there are the two solutions [5( 

ft « 0.796, #2 = 1. (18) 

A numerical solution of the set of equations without the angular approximation was pre- 
sented in ref. [3[ corresponding to the value (3\ ~ 0.796. In this paper we present an improved 
numerical method with which the infrared exponent /?2 = 1 is found as well. 

After obtaining the self-consistent solutions for the form factors numerically, we may apply 
them to the calculation of the static quark potential on the one hand, and a non-perturbative 
running coupling on the other. The static quark potential V(r) can be computed by taking 
the expectation value of the Coulomb term of the Yang-Mills Hamiltonian where p is 
chosen to be the charge distribution of two heavy color charges, i.e. 

p{x) = 5 {3) {x - r/2) - 5 {3) {x + r/2). (19) 

After Fourier transform, we arrive at 

V(r) = /|^^(l-e-). (20) 

The solution with /3 2 = 1, see Eq. (|18p . has the attractive feature that it gives rise to a 
strictly linearly rising static quark potential V(r) = or. 

A non-perturbative running coupling may be extracted from the RG invariant V(r) given in 
Eq. (|20|) and would be divergent in the infrared. At the same time, it is possible to define 
the running coupling via the ghost-gluon vertex 



16 a 2 2 
a ih) = ^.MJLk 5 G 2 (k)D(k) = —h d 2 (k) ou (k) , (21) 
3 An 3tt 

and find that, due to the non-renormalization of the ghost-gluon vertex, a(k) has an infrared 

fixed point. This can be seen immediately when plugging in the sum rule (1161) into Eq. 

(J2TJ). It was shown in ref. [5| that a(0) = for (5 = {3% by infrared analysis of the 

Dyson-Schwinger equations. The same value is actually found when resorting to the angular 

approximation as done in There, it was shown that ^ = 2 p(i3+i) w hich yields the 

above mentioned infrared fixed point if (3 = fa- 
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III. NUMERICAL SOLUTION OF THE DYSON-SCHWINGER EQUATIONS 



In order to solve the set of integral equations (jSJ) and Qj numerically, we use a Chebyshev 
approximation on a logarithmic scale for the form factors. A logarithmic momentum scale 
is introduced to have enough nodes at every order of magnitude. Furthermore, infrared- or 
ultraviolet divergent functions are expressed on a double-logarithmic scale. 
Whenever we encounter a change of sign in one of the functions, the logarithmic scale 
demands a separate treatment of positive and negative values. In particular, the function 
x{k) goes to +00 for k — > and to —00 for k — > 00, with a change of sign at some intermediate 
point which we can fix by the renormalization condition xif 1 ) = 0- This knowledge then 
allows us to fit x(h) logarithmically in the infrared, and to fit — x(fc) also logarithmically 
in the ultraviolet. For some range around k = /i, we fit linearly. The resulting set of 
nodes is to be used also for the other form factors. This is done to minimize the number of 
required Chebyshev evaluations since every evaluation of a fit is subject to some non-zero 
error. 

For asymptotic values of momentum, it is instructive for the numerics to make use of the 
analytical results, cf . the methods used in ref . [sj . The set of nodes cover a finite momentum 
region. In order to extrapolate this to the whole infinite momentum range, the general 
algebraic forms obtained analytically serve as ansatze. However, contrary to what has been 
previously done in Dyson-Schwinger studies, the parameters of these asymptotic forms are 
still determined numerically. E.g., in the infrared it was shown analytically that power laws 
provide a solution. The numerics will have the asymptotic form of the power laws ffl~5l) as 
an input but determine the various exponents and coefficients as parameters by nonlinear- 
least-squares fitting. Thus, this numerical method provides a check on the analytical results 
ffT8]) rather than imposing the latter. 

Furthermore, since the set of integral equations comprises the difference u(k) — x(fy for 
small momenta k, we extract the infrared power law A/k a of u(k) and x{k) in a single 
fit with the combined data which gives us A and a, and extract the next-to-leading order 
constants u and X01 again with least-squares fits. This gives us fits for u(k) and x(k) U P 
to intermediate momentum values, 



It is possible to gain more information on the intermediate momentum regime by setting up 
an infrared expansion of the function (cf. Eq. ( TT7j) ) 



and determining its parameters by least-squares fitting. However, this works only if the form 
factors already are pretty close to the (yet unknown) solution. Typically, at intermediate 
iteration stages u(k) is an infrared divergent quantity (as opposed to the fully iterated 
solutions). It is then mandatory to control that the fit of the function v(k) converges to a 
sensible value at each iteration step. 

Ultraviolet extrapolation is done in a similar way. We know the asymptotic behaviour from 
the UV analysis, see ref. 0] , and make corresponding ansatze for the form factors in the UV, 




(22) 



u(k) :=u(k)-x(k) = c + CiA; 7 , fc->0,7>0. 



(23) 



u{k 



d(k — > 00 



00) = k, 





(24) 



(25) 
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FIG. 1: Left panel: the gluon energy u(k) and the modulus of the scalar curvature x(k)- Right 
panel: the difference uj(k) — 



extracting the coefficients as well as the different scale parameters fj, x , fid, Hf using least- 
squares fitting. 

Equipped with the discretized form factor ansatze, we can now perform the numerical inte- 
gration by quadrature followed by iteration to determine the fit parameters. Here, we use 
the set of renormalized and thus finite integral equations © and (Q. Whereas the angular 
integrals are performed with the customary Gauss-Legendre formula, the radial integration 
involves additional non-linear transformations. Having introduced finite infrared and ultra- 
violet cutoffs e and A, resp., for the radial integrals of the kind J °° dx f(x), 1 we transform 
the integration variable logarithmically to map the integration interval to [In e, In A] . Like in 
the Chebyshev expansion, this logarithmic representation is chosen to have enough nodes at 
every order of magnitude. The remaining integral over this interval is evaluated numerically 
using a Gauss-Legendre algorithm. 

The numerical solution is then carried out iteratively, thereby using some relaxation method 
(like X new < ^X new "I - (1 r)x id) to compute the updated values for the form factors at the 
Chebyshev nodes. For a stable iteration we use r « 0.1. These new values are then used 
to extract the Chebyshev coefficients and the fit parameters for the infrared and ultraviolet 
extrapolations. This gives us new form factors. 

For the explicit calculation it is convenient to express the Dyson- Schwinger equations in 
terms of dimensionless quantities by rescaling all dimensionful quantities with appropriate 
powers of the renormalization scale /i. Indicating the rescaled quantities by bars, one has 

k = - , u(k) = U{k = ^ , d(k) = d{k = fxk) , 

a!) = ^, c=-- (26) 

Furthermore, we have to specify the renormalization constants, which are d(l), u)(l), 

c. It was already observed in ref. that neither the infrared exponents nor the ultraviolet 



1 The error thus involved can be shown to be of the order 0(e) and 0(1/ A) to some positive power. 
Choosing e sufficiently small and A sufficiently large will render the error negligible. 
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FIG. 2: Left panel: d(k). Right panel: the functions u(k) and d(k) for different values of u(fj,). 



behaviour of the various form factors depend on the precise value of the renormalization 
constants except for the value of d(l) which is, however, fixed by the horizon condition 
<i _1 (0) = 0. The remaining renormalization constants are fixed as follows. Throughout 
our calculations we choose x(l) = 0. This choice is possible since from the infrared and 
ultraviolet analysis of the Dyson- Schwinger equations carried out in ref. 0] follows that 
x(k) has a zero at some finite momentum value. Choosing x(l) = has the advantage that 
the renormalization of the curvature x(k) just removes the ultraviolet divergence without 
introducing any new finite parameter. In the same spirit we chose c = which implies that 
u(k) and x(k) approach each other for k — > (see Eq. (j!7p ). This choice is preferred by 
the study of the 't Hooft loop [l(J 11]. The remaining renormalization constant uj(1) is 
considered parameter. 

Figures [1] and [2] show the numerical solution of the coupled Dyson-Schwinger equations for 
the choice a)(l) = 1.6. As seen from the left panel in Fig. [H uj{k) and indeed obey the 
same infrared divergent behaviour, in agreement with the infrared analysis of the Dyson- 
Schwinger equations 0] . Furthermore from the right panel of Fig. [1] we also notice that 
in the infrared uj{k) and x{k) coincide, in accord with our choice of the renormalization 
condition c = 0. 

Figure [2] shows the dependency of ui{k) and d{k) upon the renormalization constant uj(h). We 
find that this constant does not change the infrared exponents and the ultraviolet behaviour 
but does have impact on the prefactors A, B of the power laws for oo(k), x{k) an d d(k). 
Also, there are some quantitative changes in the intermediate momentum range. 
Figure [3] (left panel) shows the running coupling, as defined by Eq. (I2~Tj) . It has an infrared 
fixed point which is in excellent agreement with the value predicted by the infrared analysis 
given in [1] . It is also worth mentioning that the running coupling in Fig. [3] is a monotonic 
function and has no more than one inflection point and thus does not bring up peculiarities 
in the beta function. 

In the right panel of Figure [3] we show the Coulomb potential, as given by Eq. (1201) . We find 
an exactly linearly rising behaviour (within an estimated error of less than one percent), 
again in agreement with the infrared analysis [5J|. The linearly rising potential allows us 



to fit our scale from the string tension a. Lattice calculations [12| show, however, that the 



Coulomb string tension is about a factor of 1.5 ... 3 larger than the string tension extracted 



from the Wilson loop, in agreement with the analytic result [13J that the Coulomb string 
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FIG. 3: (left) The running coupling a(k). (right) The Coulomb potential V(r). 



tension is an upper bound to the Wilson loop string tension. 



IV. SUMMARY AND CONCLUSIONS 



We have numerically solved the Dyson- Schwinger equations following from the minimization 
of the vacuum expectation value of the Yang-Mills Hamiltonian in Coulomb gauge assuming 
the horizon condition. The solution found produces a strictly linearly rising static quark 
potential which is in agreement with the analytic infrared analysis carried out in ref. [||. 
We find only a rather weak dependence of our solutions on the specific choice of the re- 
maining renormalization constants, in particular, the infrared exponents and the ultraviolet 
behaviour are independent of the choices of the renormalization constants, once the horizon 
condition is adopted. There are, however, changes in the intermediate momentum range 



which, as will be shown in ref. [lOj, are tested by the 't Hooft loop 



Acknowledgments 

It is a pleasure to thank A. Szczepaniak for valuable discussions and ongoing correspondence 
on the subject. We are grateful to the referee for pointing out ref. [9| to us. This work was 
supported by the Deutsche Forschungsgemeinschaft (DFG) under contract no. Re856/6-l,2 
and by the Europaisches Graduiertenkolleg Basel-Graz-Tiibingen (EUROGRAD). 



[1] N. H. Christ and T. D. Lee, Phys. Rev. D 22 (1980) 939 [Phys. Scripta 23 (1981) 970]. 

[2] A. P. Szczepaniak and E. S. Swanson, Phys. Rev. D 65 (2001) 025012 [arXiv:hep-ph/0107078] ; 

A. P. Szczepaniak, Phys. Rev. D 69 (2004) 074031 | arXiv:hep-ph70 306030 1 . 
[3] C. Feuchter and H. Reinhardt, Phys. Rev. D 70 (2004) 105021 [arXiv:hep-th/0408236| . 
[4] H. Reinhardt and C. Feuchter, Phys. Rev. D 71 (2005) 105002 [arXiv:hep-th/0408237| . 
[5] W. Schleifenbaum, M. Leder and H. Reinhardt, Phys. Rev. D 73 (2006) 125019 

[arXiv:hep-th/0605115] . 



S 



[6] D. Zwanziger, Nucl. Phys. B 485 (1997) 185 [arXiv:hep-th /9603203|. 

[7] J. C. Taylor, Nucl. Phys. B 33 (1971) 436; A. Cucchieri, T. Mendes and A. Mihara, JHEP 0412 
(2004) 012 | arXiv:hep-lat /0408034| ; W. Schleifenbaum, A. Maas, J. Wambach and R. Alkofer, 
Phys. Rev. D 72 (2005) 014017 larXiv:hep-ph/041 10521. 



C. S. Fischer and D. Zwanziger, Phys. Rev. D 72 (2005) 054005 arXiv:hep-ph/0504244 



[9] D. Atkinson and J. C. R. Bloch, Phys. Rev. D 58 (1998) 094036 [arXiv:hep-ph/97124 59|; 
J. C. R. Bloch, |arXiv:hep-ph/0208074 ; C. S. Fischer and R. Alkofer, Phys. Lett. B 536 
(2002) 177 |arXiv:hep-ph/0202202| ; A. Maas, Comput. Phys. Commun. 175 (2006) 167 
|arXiv:hep-ph/05041~10] . 
[10] H. Reinhardt and D. Epple, to be published. 



[11] H. Reinhardt, D. Epple and W. Schleifenbaum, arXiv:hep-th /0610324 
[12] A. 



Cucchieri and D. Zwanziger, Nucl. Phys. Proc. Suppl. 119 (2003) 727 
'arXiv:hep-lat/0209068|; J. Greensite and S. Olejnik, Phys. Rev. D 67 (2003) 094503 



arXiv:hep-lat/0302018|; K. Langfeld and L. Moyaerts, Phys. Rev. D 70 (2004) 074507 



arXiv:hep-lat/0406024 



[13] D. Zwanziger, Phys. Rev. Lett. 90 (2003) 102001 [arXiv:hep-lat/0209105| . 



9 



